#K-means算法
import numpy as np
from sklearn.datasets import load_digits
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from matplotlib.offsetbox import OffsetImage,AnnotationBbox
import matplotlib.pyplot as plt
from time import time
from sklearn import metrics # 该模块用于指标计算
from sklearn.pipeline import make_pipeline # 该模块用于多个操作连成管道统一处理
from sklearn.preprocessing import StandardScaler # 该模块用于数据归一化
from sklearn.decomposition import PCA

data, labels = load_digits(return_X_y=True)
(n_samples,n_features),n_digits = data.shape,np.unique(labels).size
print(f"#digits: {n_digits};#samples: {n_samples};#features: {n_features}")

kmeans = KMeans(n_clusters=10, n_init=4,random_state=0)
kmeans.fit(data)
y = kmeans.labels_[:100]
print(y) # 打印聚类结果前100个样本的类别标签
print(kmeans.cluster_centers_.shape) # 打印存储各类别类中心的数组的形状

pca = PCA(n_components=2) #实例化，定义一个PCA模型，选择前两个主成分
X_r = pca.fit_transform(data[:100]) # 值对前100个样本进行降维用于可视化

fig, ax = plt.subplots(1,2,figsize=(10,4)) # 分左右两个子图显示可视化效果
# 设置左子图绘制颜色
colors = ['navy','turquoise','darkorange','red','blue','yellow','black','green','grey','#123456']
# 设置右子图绘制颜色，因左右子图绘制方式不同，颜色不能通用
cmapcolors = ['Greens','Greys','pink_r','Purples','gist_earth_r','PuRd','Oranges','PuBu','gray_r','BuGn']
# 绘制左子图，用不同颜色的点代表样本
for color, i in zip(colors,range(10)):
      ax[0].scatter(X_r[y==i,0],X_r[y==i,1],color=color,alpha=0.8,lw=2,label=str(i))
      ax[0].legend(loc="best",shadow=False,scatterpoints=1,fontsize='xx-small')
      # 绘制右子图，直接显示不同颜色的数字图像
ax[1].set_xlim(-35,35)
ax[1].set_ylim(-30,30)
for i in range(100):
      imagebox=OffsetImage(data[i].reshape(8,8),zoom=1,cmap=cmapcolors[labels[i]],alpha=0.65)
      ab = AnnotationBbox(imagebox,xy=X_r[i],frameon=False)
      ax[1].add_artist(ab)
# 可视化
plt.show()

# 通过参数init来设置三种不同的聚类中心初始化方法，对KMeans算法进行性能测试。
print(69 * "_") # 打印分隔线
print("init\t\ttime\tinertia\thomo\tcompl\tv-meas\tARI\tAMI") # 打印表头
print(69 * "_")
pca = PCA(n_components=10).fit(data) # 样本降维到10
# 以下分别定义5种聚类初始方法的参数
init_name = ['k-means++(1)','k-means++(4)','random(1)','random(4)','PCA']
init_param = ['k-means++','k-means++','random','random',pca.components_]
n_init_v = [1,4,1,4,1]

for n_init, name, init in zip(n_init_v,init_name,init_param):
      # 应用不同参数进行实例化
      kmeans = KMeans(init=init, n_clusters=n_digits,n_init=n_init,random_state=2)
      # 计算此时模型各项指标并输出
      t0 = time()
      # 使用管道执行聚类算法
      estimator = make_pipeline(StandardScaler(),kmeans).fit(data)
      fit_time = time() - t0 # 得到训练时间
      results = [name,fit_time,estimator[-1].inertia_] # 第3项为类内总距离平方和
      # 定义需要计算的性能指标，以下指标值越大代表性能越好
      clustering_metrics = [metrics.homogeneity_score, # 同质性度量
                            metrics.completeness_score, # 完整性度量
                            metrics.v_measure_score, # 调和平均值
                            metrics.adjusted_rand_score, # ARI指数
                            metrics.adjusted_mutual_info_score #AMI指数
]
      # 利用原始标签和预测得到的标签得到每个性能指标值
      results += [m(labels,estimator[-1].labels_) for m in clustering_metrics]
      # 展示评估指标结果
      formatter_result = ("{:9s}\t{:.3f}\t{:.0f}\t{:.3f}\t{:.3f}\t{:.3f}\t{:.3f}\t{:.3f}")
      print(formatter_result.format(*results))
print(69 * "_")